function [J grad] = linearRegCostFunction(theta, X, y, lambda)

    m = size(X, 1);
    X = [ones(m, 1), X];
    h = X * theta;
    theta_reg = [0; theta(2:end)];
    J = 1 / (2 * m) * sum((h - y).^2) + lambda / (2 * m) * theta_reg' * theta_reg;
    grad = X' * (h - y) / m + lambda / m * theta_reg;

end
